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In a bilayered system of particles with wake-mediated interactions, the action-reaction symmetry 
for the effective forces between particles of different layers is broken. Under quite general conditions 
we show that, if the interaction nonreciprocity exceeds a certain threshold, this creates an active 
dispersion of self-propelled clusters of Brownian particles. The emerging activity promotes unusual 
melting scenarios and an enormous diffusivity in the dense fluid. Our results are obtained by com¬ 
puter simulation and analytical theory, and can be verified in experiments with colloidal dispersions 
and complex plasmas. 
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I. INTRODUCTION 

Effective forces between mesoscopic particles often be¬ 
come nonreciprocal when the interactions are mediated 
by a nonequilibrium environment. Such situations can 
be realized in various soft matter systems - most not ably 
in colloidal dispersions [H-Q and complex plasmas (^10l| . 
where microparticles are embedded, respectively, in a liq¬ 
uid solvent or a dilute weakly ionized gas. In particular, 
the action-reaction symmetry in these systems is broken 
when the surrounding fluid moves with respect to the 
particles (lll - [l5l| . or when the interaction of molecules 
with the particle surface is out of equilibrium Qiii. 

Studies of nonreciprocal interactions have gained in¬ 
creased interest in recent time. When the dynamics of 
individual particles is undamped (Newtonian) or weakly 
damped (when the relevant dynamical timescales are 
much shorter than the damping time) (l^ . which is typ¬ 
ical for complex plasmas, one can observe a remark¬ 
able state of detailed dynamic equilibrium with different 
species having different temperatures. For Brownian dy¬ 
namics, it has recently been shown that mixtures of dif- 
fusiophoretic colloids experience effective nonreciprocal 
forces which stimulate the formation of stable aggregates 
(so-called active molecules) and trigger collective os¬ 
cillatory motion [l^ . 

In this paper, we consider a broad and generic class 
of nonreciprocal interparticle forces, the so-called wake- 
mediated interactions. Particles embedded in a flowing 
medium generate wakes, which contribute to the inter¬ 
actions with neighbors in a nonreciprocal way im. 
Similarly, particles emitting chemicals in a certain di¬ 
rection generate asymmetric concentration fields - artifi¬ 
cial “chemical wakes”, also exerting non-reciprocal forces 
on the neighbors [f^ 113 • In both cases, the action- 
reaction symmetry is only restored for identical particles. 
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forming a perfect monolayer perpendicular to the direc¬ 
tion of wakes. Therefore, here we consider a quasi two- 
dimensional system of Brownian particles [f8l | which are 
kept into stable bilayers by external fields such as elec¬ 
tric, magnetic, gravitational or optical fields. Under quite 
general conditions posed on the mutual reciprocal and 
nonreciprocal forces, we observe a continuous transition 
from inactive (stacked) pairs to active units, indicating 
the emergence of active fluids. Different from ordinary 
active particle systems , these active units can 

break and become passive again. Using analytical the¬ 
ory and simulation including hydrodynamic interactions 
between the particles, we explore the full density regime 
up to freezing and find an unusual melting upon densifi- 
cation, along with a reentrant freezing and an enormous 
diffusivity in the concentrated fluid. 

The paper is organized as follows: In Sec. [TTlwe specify 
our model, perform a stability analysis for small clusters 
in Sec. lIIIl and describe our simulation in Sec. lIVI Results 
are discussed in Sec. El and summarized in Sec. EB 

II. MODEL 

The motion of a particle i at position in the plane 
is governed by the fully damped Langevin equation 

f. = y:L„ (F,+{,) + itBry:^, (1) 

j 3 ^ 

where is the mobility matrix and is a random force. 
The total force on particle i is given by F^ j, 

where F^j is the pair-interaction force exerted by a par¬ 
ticle j on the particle i. The random force is Gaussian 
distributed with zero mean, (^^(t)) = 0, and variance 
{$,i{t)$,j{t')) = 2 k-QT5{t — t'), where T is the ther¬ 
mostat temperature, /cb the Boltzmann constant, 5{t) the 
Dirac delta function, and the inverse of L. In this 
paper, we include hydrodynamic interactions in the zero- 
temperature limit, and neglect them at finite tempera- 
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tures. The latter approach is justified when the suspen¬ 
sion is highly dilute, but still strongly interacting. Then, 
each mobility matrix reduces to Ly = 7 ~^( 5 yT, with the 
unit matrix I and a friction coefficient 7 ^. In the zero- 
temperature limit, we consider the mobility matrix to be 
approximated by the Oseen tensor [2lj 

^oir) = ^il + rr), ( 2 ) 

where Rh is the hydrodynamic radius, r = |r| and f = 
r/r. Thus, we have Ly Ri Lo(ri — r^) for i ^ j and 
Liz — 1/Ti ■ 

We consider a typical situation when interactions be¬ 
tween particles are isotropic in the plane. In this case 
the mutual forces between particles i and j are radial, 
i.e., Fy = FijUiij with riy being the unit vector from 
i to j, and Fy only depends on the absolute distance 
Ty = \ri — rj\. Furthermore, we introduce species A 
and B and attribute particles to the same species if their 
pair interactions are reciprocal, i.e., FAA('f) = FBB('r) = 
—dipr{r)/dr. A generic form for the forces between dif¬ 
ferent species, 

FAB.BA(r') = -dipr{r)/dr ± d(pn(r)/dr, (3) 

is a superposition of the reciprocal (r) and nonreciprocal 
(n) components, determined by the respective potentials 
(^r,n- The latter are related to the potential t/jy generated 
by the particle i at the location of the particle j via ipr,n = 
± ‘Pij)- Thus, the pair interactions are reciprocal if 
ifij = ifji, and are nonreciprocal otherwise. 

An important class of a constant nonreciprocity is re¬ 
alized when ipr{r) and (pnir) are similar functions, i.e., 
when the nonreciprocity = A = const. For 

the undamped (Newtonian) or weakly damped dynamics 
with A = const, the equations of motion can be equiv¬ 
alently transformed into a reciprocal form by a simul¬ 
taneous proper renormalization of the interaction forces 
and masses, i.e., such dynamics can in fact always be 
described by a (pseudo) Hamiltonian [l^. It is note¬ 
worthy that for the Brownian dynamics with nonrecipro¬ 
cal interactions one can employ a similar approach: By 
renormalizing the interactions with Eq. (3) of Ref. (l^ . 
and introducing the renormalized damping coefficients 
7a. B = 7a.b/(1 T A), we readily transform Eq. ([T]) to 
the form where the interactions are reciprocal, while the 
solvent temperatures for different species A and B are dif¬ 
ferent and equal to Ta.b = T/( 1 =f A). Interestingly, such 
a “hetero-Brownian” model has been recently introduced 
in a different context, to describe DNA dynamics 
and was also proposed for colloidal pairs under external 
forcing (23 - [^ . 

In our model, we adopt a binary mixture of point¬ 
like particles whose direct (electrostatic) interactions are 
characterized by the charges Qa and Qb- The particles 
are confined in a horizontal xy-plane in two layers with a 
height difference h, as sketched in Fig. [T] The point¬ 
like approximation is justified as long as the distance 




FIG. 1. Schematic sketch of nonreciprocal wake-mediated in¬ 
teractions. The particle species A and B are confined in the 
upper and lower layers, respectively. While the direct inter¬ 
particle forces are reciprocal, -|- F^g = 0, the particle- 
wake forces are nonreciprocal, Fg^ + F^b 7^ 0 such that for 
the total forces Fab 7^ —Fba. 


between the particles is much larger than their diame¬ 
ter. For non-reciprocal particle-wake forces we have two 
representative realizations in mind, both leading to effec¬ 
tive Yukawa interactions: The first one can be obtained 
with catalytically active Janus colloids (^ . [ 2 ^, where 
each particle emits a chemical from the lower segment of 
the surface, thus generating “chemical” wakes (the verti¬ 
cal orientation of such particles can be stabilized in an 
external magnetic field [^). The emitted chemicals de¬ 
compose with a constant rate , such that the effective 
interaction between the chemical wake and a neighboring 
particle is of Yukawa form. The second realization is an 
externally imposed micro-ion flow, parallel to the vertical 
2 -axis, which induces an “electrostatic” wake below each 
particle, while the fluid remains at rest. In both cases, 
the wake is mimicked by a point-like effective “charge” qi 
at the distance S below each particle. 

The total interaction force is a combination of the 
direct particle interaction and the particle-wake inter¬ 
action. The force depends on several free parameters 
whose combination determines emergence of different 
self-organization phenomena. In this Letter we chose a 
certain set of parameters, demonstrating variety of the 
emerging activity, and identify a general necessary con¬ 
dition for the activity onset [see Eq. (O below]. Let 
us introduce the three-dimensional particle coordinates 
Ri and the corresponding coordinates in the horizon¬ 
tal plane. Then the force F^ = —dipijjdrj exerted in 
the horizontal plane by the particle i on the particle 
j is determined by the potential tpij = QiQjY[Rij) + 
diQjY{RJj), where Y{R) = is the (unity 

charge) Yukawa potential which depends on the distance 
Rij = |Ri — Rjj between the particles as well as on the 
distance = |Rj —R^ —between the particle j and 
the wake center of the particle *; here, we assume that a 
Yukawa potential for both forces has the same effective 
screening length A, and that qi oc —Qi [3ll - [3^ . For parti¬ 
cles in the same layer, A or B, we have i?7 = there- 
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fore, (/9n = 0 (since qiQj = qjQi), and hence the forces 
are reciprocal. For the AB interactions the symmetry is 
broken, ^ RJ^, and the forces are nonreciprocal. For 
simplicity, particles A and B have charges of the same 
magnitude and opposite signs, Qa = —Qb = Q, the 
same friction coefficients, 7 a = 7b = 7) and the height 
difference is ft, = A. A natural measure of nonreciprocity 
in this case is the relative wake charge, q = —qijQi > 0. 


III. STABILITY ANALYSIS FOR SMALL 
CLUSTERS 

In order to illustrate a tendency of particles with non¬ 
reciprocal interactions to self-organize themselves with 
increasing q, and to identify the characteristic building 
blocks of this complex process, let us consider the for¬ 
mation of small clusters in the absence of hydrodynamic 
interactions. Then, the equilibrium configurations for a 
cluster of N particles are determined from the force bal¬ 
ance in the horizontal plane, 

N 

^F,,(r,,)=F, (4) 
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where the net force F is a constant horizontal vec¬ 
tor for Vi € We apply the standard stabil¬ 

ity analysis of the derived configurations in the zero- 
temperature limit. This corresponds to the eigenvalue 
problem det (9Fy|eq — ywl) = 0, where ... |eq de¬ 
notes the {2N X 2N) dynamical matrix calculated for the 
equilibrium configurations. 


A. Doublets 

A pair of particles of different species form an equi¬ 
librium doublet with the horizontal separation ro when 
Fab(i'£i) = —Fba(i'£>) = F; the doublet is stable if 
d [Fab(i’) + Fba(i’)] ® 

elude that the stability condition is only fulfilled when 
the reciprocal component of the force is equal to zero, 
d(fr{r)/dr\^^^^ = 0 . 

For a vertical pair ro = 0 - we call it an inactive 
doublet - two regimes can be distinguished: (i) When 
dFij{r)/dr\^^g < 0 for both particles, they return to the 
equilibrium after a small perturbation. Below we demon¬ 
strate that this case, sketched in Fig.[2](a), is observed for 
a “weak” nonreciprocity, when the relative wake charge 
is smaller than a certain critical value, q < gcri (i-e-, 
this always occurs for reciprocal interactions), (ii) When 
dFij{r)/dr\^^g > 0 for one of the particles, the restor¬ 
ing forces are pointed in the same direction, as shown 
in Fig. HKb). The equilibrium in this case, correspond¬ 
ing to dcri < q < qcT 2 , would only be restored in the 
zero-temperature limit; in the presence of an infinitesi¬ 
mal thermal noise the doublet should break apart. 


FIG. 2. Stable configurations of two-particle clusters, depend¬ 
ing on the relative wake charge q. (a-c) Sketches illustrate 
three distinct regimes (side view): For q < gcri, particles form 
a stable vertical pair, an inactive doublet, since the restoring 
forces Fab and Fba exerted by a small perturbation pull the 
particles back; for qai < q < qci 2 , the vertical pair remains 
stable only in the zero-temperature limit assumed here, since 
Fab and Fba are pointed in the same direction; for q > qci2, 
the particles form an active doublet with a finite horizontal 
separation, moving along the force Fab = Fba. (d,e) Equi¬ 
librium horizontal separation of the doublet, ro (normalized 
by A) and the corresponding doublet velocity vd (normalized 
by Q^/A^ 7 ), the shading indicates the stability regimes illus¬ 
trated in (a-c). The results are for the wake length <5 = 0.2A, 
ft = A and Rh = 0 (black dashed line) as well as Rh = 0.2A 
(blue solid line). 

Under the general condition 

diprir)/dr\^^^^ =0 and dipnir)/dr\^^^^ ^0, (5) 

satisfied for rqj > 0, a pair emerges which is self- 
propelled in the direction hab with the velocity vd = 
— 7 “^ dipn{r )/Such clusters will be referred to 
as active doublets and occur when q > gcr 2 - Note that for 
a constant nonreciprocity, A = const, stable doublets are 
always at rest, since (pn{i") = A(^r(r) and therefore the 
nonreciprocal force is equal to zero at r = r^i . Equa¬ 
tion ([S]) represents the necessary condition for the emerg¬ 
ing activity, which can be satisfied for various combina¬ 
tions of the interaction parameters (e.g., when Qa and 
Qb have the same sign, but the direct and the particle- 
wake interactions are characterized by different screening 
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lengths). 

Figures [DJd) and (e) illustrate the results of the sta¬ 
bility analysis in the horizontal plane, performed in the 
zero-temperature limit. In the present example, two par¬ 
ticles of different species are stacked on top of each other 
(i.e., they form an inactive doublet) when the relative 
wake charge is smaller than gcr2 — 0 . 74 . For larger q, 
the separation rjy continuously increases and an active 
doublet moves along its symmetry axis, with the velocity 
vjy which varies non-monotonically with q. Thus, in di¬ 
lute systems (with infinitesimal number density) one can 
expect the formation of multiple individual doublets. Hy¬ 
drodynamic interactions do not influence the pair separa¬ 
tion ru, but cause a velocity increase. Interestingly, the 
leading term for the hydrodynamic far-held for an active 
doublet is a force monopole, as opposed to a standard 
microswimmer where it is a force dipole (^ . 

For a hnite number density p (number of particles per 
unit area) in systems without hydrodynamic interactions, 
we analyze the stability of crystalline structures in the 
zero-temperature limit. The time-dependent coordinate 
of the ith-particle is presented as a sum of its equilibrium 
lattice position and a displacement, ri{t) = r^q^i + 

The interaction force, Eq. ®, is then expanded to the 
hrst order in and substituted in Eq. ©• Using 
Ui oc exp(ik • req,i + cut), the dispersion relations u;(k) 
are derived as eigenvalues of the resulting dynamical ma¬ 
trix ■ We examine a vertically stacked hexagonal lat¬ 
tice and an interdigitated hexagonal lattice, the stability 
requires Re w(k) < 0 for all k from the hrst Brillouin 
zone of the lattice. 

Let us now study the effect of hydrodynamic interac¬ 
tions on a doublet, consisting of a particle of species A 
and one of species B in the dilute limit. The general 
requirement for stable clusters is the equality of the ve¬ 
locities, 


rA = re ■ 


( 6 ) 


The equation of motion of the A particle can be written 
as 


i-A = -FA(rA - re) + 
7 


3 Rg 

dyf 



FB(rB - va ) , 


( 7 ) 

with r = + h'^', the respective equation for the B 

particle, is obtained by the A •«->• B permutation. Equa¬ 
tion dni) is only fulhlled if Fab(?'d) = — Fba(?'d), as in the 
case without hydrodynamic interactions. Thus, Eqs. m 
remain valid also for a hnite hydrodynamic radius. Fig¬ 
ure [ 5 ] shows the result for a doublet with and without 
hydrodynamic interactions. The doublet distance rjj re¬ 
mains unchanged, while the doublet velocity is slightly 
increased. 


B. Triplets 


inactive singlet inactive singlet and 

and doublet active doublet 



active linear inactive linear active triangular 
triplet triplet triplet 


FIG. 3. Stable configurations and velocities of three-particle 
clusters in the zero-temperature limit. The two possible com¬ 
binations of the species A and B are shown. The figure legend 
is the same as in Fig. [2] 


composed of one particle B and two particles A with 
negligible hydrodynamic interactions. We work in the 
frame of reference of the hrst (B) particle, i.e., the coor¬ 
dinates of the second and third (A) particles are r2,3. In 
this case, the general equilibrium condition, Eq. 0 , can 
be identically transformed to the following two equations 
for the particle coordinates (plus one equation for the net 
force F): 

F(r2)n2 -I- F(r3)n3 = 0 , ( 8 a) 

2 FAA(r 23 )n 32 — FBA(r3)n3 + FBA(r2)n2 = 0, (8b) 

with F{r) = Fba('c)+ 2 Fab(?')- In the reverse case, where 
clusters are composed of one A and two B particles, the 
labels are simply to be swapped. 

Using Eq. ((8al) . one can distinguish two principal cases: 

(i) F(r2 3) ^ 0, then solutions exist only for r2 || r3; 

(ii) F{r2) = F{r^) = 0, then solutions are possible for 
noncollinear r2 and r3. 

(i) If F{r) is a monotonic function, the only solution 
is r2 = — r3; from Eq. ()8bl) we obtain FAA( 2 r) = 
FBA(r), which yields r2 = = r. Due to symme¬ 

try, F = 0 and hence we call such conhgurations 
inactive linear triplets. However, if F{r) is a non¬ 
monotonic function, also solutions with r2 7^ r^ 
are possible - in this case Eqs. (I 5 al) and (I8bl) are 
reduced to F{r2) = Firs) and 2 Faa(?’2 + rs) = 
Fba(?'2) + Fba(?' 3)- Such asymmetric clusters usu¬ 
ally imply a non-vanishing net force, F 7^ 0, which 
generates a directed propulsion. We call these con¬ 
figurations active linear triplets. 


Eor three particles, there is a variety of possible triplet 
configurations. To start with, let us consider a cluster 


(ii) If Fba(?') is monotonic, solutions for r2_3 are lim¬ 
ited to triangles with r2 = r^ = r and the 
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Na = 2 
Nb = 1 




r2 



FIG. 4. Equilibrium horizontal separation in the doublet and 
triplet configurations shown in Fig. [3l 


9 < 180°. 

If hydrodynamic interactions are taken into account, 
the equilibrium condition generalizes towards 


Ti =r V i G [1, TV] . (9) 


Generally, the triplet coordinates derived above do not 
fulfill this equilibrium condition, which is in contrast with 
the doublets (where the inclusion of hydrodynamic inter¬ 
actions only induces a rescaling of the velocity). For three 
particles, we solve Eq. m numerically and show the re¬ 
sults in Figs. O and m However, the resulting changes to 
the particle coordinates are minor, as shown in Fig. |4l 
Similar to doublets, the inclusion of hydrodynamic inter¬ 
actions merely causes a slight increase of the velocity, see 
again Fig. [31 


apex angle 9, obtained from F{r) = 0 and 
FAA(2rsini0)/sini0 = Fba(»’)- Such configura¬ 
tions are called active triangular triplets. Finally, 
if Fba(^) is a non-monotonic function, triangular 
triplets with r 2 ^ rs are possible. 

As for the doublets, we apply the standard stability anal¬ 
ysis of the derived configurations in the zero-temperature 
limit. 

If the number of particles B is twice as high as the num¬ 
ber of A particles, the dependence on q remains the same 
as in Fig. |2l The “excess” particle B simply remains an 
inactive singlet [see Fig. |3][a)]. On the contrary, in the 
situation with two particles A for one particle B vari¬ 
ous active and inactive structures emerge, as presented 
in Fig. |3Kb): An inactive doublet and an inactive sin¬ 
glet are formed when q < 0.17, while for q G (0.17, 0.36) 
they merge into an active linear triplet, where the po¬ 
sition of particle B is slightly shifted from the center 
(which determines the propagation direction along the 
symmetry axis). The linear triplet becomes inactive at 
q G (0.36, 0.74). The further increase of the wake charge, 
q G (0.74,0.81), causes particle B to shift perpendicu¬ 
lar to the symmetry axis, leading to an active trian¬ 
gular triplet. For even larger values of q, the triplet 
breaks apart and an active doublet and an inactive singlet 
emerge. In a similar manner, one can straightforwardly 
generalize the analysis for larger clusters or investigate, 
e.g., the rotation activity. 

Figure Ufa) shows the horizontal separation ri for the 
case of {Na = 1,Nb = 2 ), where a passive singlet and a 
doublet form. In the reverse situation {Na = 2, Nb = 1), 
we characterize the emerging triplets by their individ¬ 
ual bond distances r 2 ,r 3 [Fig. |D(b)[ and the respective 
apex angle 9 [Fig. |3Kc)[. Activity is a result of symmetry 
breaking, therefore active units are found if r 2 ^ ^3 or 


IV. COMPUTER SIMULATIONS 

We solve the equation of motion, Eq. o, using a 
forward time-step algorithm in a Brownian dynamics 
simulation for three distinct cases: when the hydro- 
dynamic interactions are neglected, we consider (i) the 
zero-temperature limit and (ii) finite temperatures; (iii) 
the effect of the hydrodynamic interactions is studied in 
the zero-temperature limit. We use a 2D rectangular 
simulation box with periodic boundary conditions and 
the edge ratio Ly/L^ = The particles are ini¬ 

tialized on a distorted stacked hexagonal lattice with a 
fixed number density p = N/{L^Ly). In the case (i) 
and (ii), we use N = 2 x 2500 particles. The respec¬ 
tive edge lengths of the simulation domain are vary¬ 
ing from {L,c,Ly) ~ (240A, 210A) at low densities to 
{Lx,Ly) ~ (43A, 38A) at pX^ = 3. In the case (iii), we 
use iV = 2 X 576 particles, and the simulation domain 
is between {Lx,Ly) ~ (115A,100A) at low densities and 
{Lx,Ly) ~ (28A, 24A) at high densities. The modeling 
of the hydrodynamic interactions is done on the Oseen 
level with Rh = 0.2A, i.e. up to volume fractions of about 
15%. For all cases, the wake length is d = 0.2A, the height 
between the layers is h = A, the time t is measured in 
units of r = qA^/Q^ and the distance r in units of A. We 
set the time step to 6t = 0.005r in the cases (i) and (ii) 
and to 6t = 0.0025r in case (iii), which ensures proper 
resolution of the particle dynamics. After initialization, 
the system is given time of IO^t to relax into a steady 
state. Statistics is gathered for multiple simulations runs 
with independent initializations and the simulation time 
of 2500t. By measuring the displacement of individual 
particles within the time step, a particle velocity is cal¬ 
culated as Vi(t) = [ri{t -b St) — ri{t)]/6t. 
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FIG. 5. State diagram in the zero-temperature limit, plotted 
in the plane of the number density p and relative wake charge 
q. Color coding depicts results obtained from the stability 
analysis, symbols show numerical results. Inactive systems 
(-I-) can be either stacked hexagonal solid (green background) 
or interdigitated hexagonal solid (blue background). For ac¬ 
tive fluid regimes (Q, red background), the average particle 
velocities are indicated by a gray scale. Diamonds (<0>) are 
used instead of circles if active doublets emerge whose decay 
time td exceeds a threshold of lO^r. The states are illustrated 
by typical snapshots, see also movies in the Supplemental Ma¬ 
terial [ 35 II . 


V. RESULTS 


A. Zero-temperature limit 

The above analytical results, see Sec. Ell are comple¬ 
mented with a numerical analysis Figure [5] presents 

the state diagram of the emerging activity, where we 
compare the theoretical results against simulations in the 
zero-temperature limit. The state diagram is plotted in 
the plane spanned by number density p and relative wake 
charge q. We identify three distinct domains: Toward 
the reciprocal limit q = 0, the particles form a bilay- 
ered stacked hexagonal crystal (green); for larger q, at 
increased density the system goes into an interdigitated 
hexagonal solid (blue); for even larger q, at low densities 
we find an active regime where the crystal melts (red). In 
addition, we show the results obtained from the numer¬ 
ical simulations. Here, we differentiate between inactive 
solids (-b) and active fluids (O)- The mobile units of 
the fluid are active doublets, that behave similar to (de¬ 
formable) active Brownian particles [stI - I^ . The active 
regime in the simulations is defined for the average parti¬ 
cle velocity (|v|) above a threshold of 10“^A/r. One can 
see that the emerging state diagram exhibits a reentrant 
behavior both with p and q. Notably, for intermediate 9 , 
there is an anomalous “water-like” melting upon an in¬ 
crease in p followed by reentrant freezing. Furthermore, 


= 0.10 -A-pA^ = 0.50 -B-pA^ = 1.00 -▼-pA^ = 1.75 
i?H/A = 0.0 i?if/A = 0.2 



0.5 0.6 0.7 0.8 0.9 1.00.5 0.6 0.7 0.8 0.9 1.0 

q q 


FIG. 6. Gharacteristics of active fluids: (a) and (d) average 
particle velocity (|v|), (b) and (e) alignment parameter c, and 
(c) and (f) decay time of doublets td, plotted versus the rela¬ 
tive wake charge q for several values of the number density p. 
The panels on the left show the results for simulations with¬ 
out hydrodynamic interactions, on the right hydrodynamic 
interactions are included. The dashed line in panel (a) shows 
the velocity of a single active doublet vd in the dilute case 
(pA^ 1), the horizontal lines in panels (a), (c) and (f) indi¬ 
cate the threshold values of the velocity (10“^A/r) and decay 
time (lO^r). 


to quantify the stability of doublets we define N]o{t), 
the average number of particle pairs that remain near¬ 
est neighbors over the time interval t. Generally, it is 
well described by an exponential decay, Njoit) oc 
with a doublet decay time T]j. Long-living active clus¬ 
ters are marked by a diamond in Fig.E The existence of 
a finite decay time td reveals a qualitative difference of 
our system to a system of permanently active particles 

[E[2flli0,li3. 

We now discuss the characteristics of the active fluid in 
more detail. Let us first introduce the averaged velocity 
(v) = {[vi{t-\-5t) — ri{t)\/5t), seeFig.jHKa) and (d), as well 
as an alignment parameter c = |(v)| / (|v|), see Fig.jSDb) 
and (e): c = 1 for a perfect nematic order and c = 0 
in a totally disordered case. The stability of doublets 
is quantified by No{t), the averaged number of particle 
pairs that remain nearest neighbors over the time interval 
t. The doublet decay time tjo is shown in Fig. [SKc) and 
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— g = 0.8 —q = 0.9 —q = 1.0 



FIG. 7. Angle of reflection Or versus the incidence angle Qi for 
a collision of two active doublets, plotted for different values 
of the relative wake charge q. The symmetric scattering, ar = 
ai, is indicated by the dotted line. 


(f). If no doublet splits during the simulation time of 
2500t, then td is set to infinity. 

Figure [B]Ja) demonstrates that at low densities (pA^ = 
0 . 1 ), the average velocity (|v|) is well reproduced by the 
velocity of a single active doublet, vd, calculated ana¬ 
lytically. Above the threshold value of gcr 2 = 0.74, the 
distance rn increases, see Fig. [2])d). For this reason, the 
average velocity hrst increases with q, but then it starts 
falling off due to decreasing interaction strength of a dou¬ 
blet, see also Fig. He). As the activity sets in, long-living 
doublets are formed throughout the system and their mu¬ 
tual collisions lead to the velocity alignment, see Fig. H 
since the angle of reflection after their mutual collision 
is always smaller the the incidence angle a^, as shown in 
Fig. [71 With increasing the number density p the onset 
of activity shifts towards smaller q, whereby the aver¬ 
age velocity vs. wake charge becomes a non-monotonic 
function, leading to a reentrant effect for pX^ > 1.25, 
where an inactive interdigitated hexagonal solid emerges, 
see Fig. [S] 

The effects of hydrodynamic interactions are demon¬ 
strated in the right panel of Fig. [51 One can see that the 
effects become significant at higher densities. In com¬ 
parison with the left panel, we observe a drastic velocity 
increase, while the alignment remains strong at any q. 
These effects become more pronounced due to the long- 
range nature of hydrodynamic interactions, so that more 
particles are involved in the collective motion. 


♦pA^ = 0.01 ■■■pA^ = 1.00 -^pA^ = 1.75 



0.0 0.2 0.4 0.6 0.8 1.0 

q 


FIG. 8 . Ratio of the long-time to short-time diffusion coef¬ 
ficients, Dl/Ds, obtained from the time dependent diffusion 
coefficient D{t) for a finite temperature T = 10“® /{k^X). 

The shading indicates a transition between active fluids 
{Dl/Ds > 1) and solids {Dl/Ds < 1). The inset depicts 
the normalized D{t) for pA^ = 0.01 and q G {0, 0.5, 0.6, 0.8}, 
demonstrating the transition from subdiffusive to ballistic in¬ 
termediate behavior with increasing q. The dashed lines rep¬ 
resent the analytical solutions of the Langevin equation for 
the diffusion of a doublet in the dilute limit, Eq. (1111) . the 
solid lines and symbols show numerical results. 


B. Finite temperature study 

Let us study the case of finite temperatures in the ab¬ 
sence of hydrodynamic interactions. For a single stacked 
doublet, where = 0, we compute the mean-squared 
displacement from Eq. (HD- We take the Taylor expansion 
of the forces around the equilibrium positions. Then, the 
radial force perturbations are SFBA{t) ~ Ca [i5rA(t) — 
(5rB(t)] and SFAB{t) Cb [SrB{t) - SrA{t)], where Ca 
and Cb are the prefactors of the linear-order terms in 
the expansion. We define A as a matrix containing these 
prefactors, such that Eq. © can be written as 

7^=A(g)X(t)+T(t), (10) 

where X(<) = {(5rA(t), (5rB(t)} is the vector of the par¬ 
ticle positions and T(t) = {fi(f), ^2(0} the vector with 
the random (radial) forces acting on the particles. For 
simplicity, we set the friction coefficient 7 independent 
of the particle index. Using variation of constants, this 
differential equation is solved by the integration over a 
matrix exponential: 

X(t) = - f dr exp [A (t - t)/ 7 ] T(r), 

7 Jo 

with (T(<)) = 0 and {T,{t)Tj{t')) = 2'ykBT5ijd{t - t')l, 
leading to (X(t)) = 0. By computing the mean squared 














-X- 0.001 


0.60 


keT A/Q2 



q 



FIG. 9. (a) Ratio of the long-time to short-time diffusion 

coefficients, Dl/Ds, obtained from the time dependent dif¬ 
fusion coefficient D{t) for given finite temperatures and num¬ 
ber density pX^ = 0.01. The horizontal dotted line indicates 
the transition between active fluids {Dl/Ds >1) and solids 
{Dl/Ds < 1). The dashed line represents the analytical solu¬ 
tion, the solid lines and symbols show numerical results, (b) 
Temporal evolution of D{t) normalized by Ds for two chosen 
relative wake charges q = 0.2 and q = 0.7 with pX^ = 0.01 
and varied finite temperature. 


displacement, we determine the diffusion ratio, 

DL{q) Cl{q)+CUq) 

Ds [CA{q) + CB{q)f' 

The result of Eq. (fTT|) is presented in Figs. |5] and IHl by 
the black dashed line. 




FIG. 10. (a) Temperature dependence of the onset of activity 
(jcri (determined from the condition Dl/Ds > 1) in simu¬ 
lations with pX^ — 0.01. At low temperatures it approaches 
the analytically derived value of (jcri — 0.46. (b) Temperature 
dependence of the doublet decay time td for simulations with 
pX^ — 0.01 and q — 0.7. If no pair splits during the simula¬ 
tion time, Td is set to infinity. The horizontal line shows the 
stability criterion of long-living active clusters. 


Finally, we demonstrate the impact of finite temper¬ 
atures in many-body simulations. In Fig. |S1 we com¬ 
pare the analytical results for the diffusion coefficient, 
Eq. HU), with the numerical results obtained for number 
density pX^ = 0.01 using the time-dependent diffusion 
coefficient D{t) = ^(|r(t) —r(O)p) plotted in Fig.[8]satu- 
rates towards a long-time diffusion coefficient Dl at long 
times which is naturally normalized by the short-time 

coefficient Ds = lim Il(t) = fc^T/y. 

(->■0 

For intermediate times there is either a sub-diffusive 
regime due to particle caging, or a ballistic regime aris¬ 
ing from the emerging activity (23, As discussed 

above (see Fig. [2]), in dilute systems the activity at finite 
temperatures is expected to set in at g > gcri(— 0.46). 
From Fig. [8] we see that for pX^ = 0.01 the transition to 
active fluids, Dl/Ds > 1, indeed occurs near this value. 
The long-time diffusion increases over several orders of 
magnitude as a function of nonreciprocity q. Even in 
dense colloidal fluids (at pX^ = 1.75) the ratio Dl/Ds 
exceeds 5, implying that there is an enormous diffusivity 
relative to the case of infinite dilution. As revealed by 
the snapshots in Fig. [SI this is mainly due to significant 
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local alignment in the fluid, which allows for an efficient 
traveling of active doublets. 

One can see that the onset of activity, defined by 
the diffusion ratio D^/Ds > 1, remains practically un¬ 
changed at all temperatures, whereas the asymptotic de¬ 
viation of D{t)/Ds from unity decreases with T, see 
Fig. O Naturally, at higher temperatures the stability 
is gradually decreasing. The explicit dependence of the 
activity onset {(jcri) and pair stability (td) on the tem¬ 
perature is shown in Fig. 1101 

VI. CONCLUSION 

In conclusion, we have shown that in two-dimensional 
systems with wake-mediated interactions a rich vari¬ 
ety of self-organization phenomena occur. In the zero- 
temperature limit. The nonreciprocal forces exerted by 
wakes generate a complex diagram of steady states. In 
particular, we showed the formation of active units - 
bound particle pairs, having interesting similarities with 
permanently active Brownian particles - and the realiza¬ 


tion of unusual melting scenarios. At finite temperatures 
we identified regimes of anomalously high diffusion. The 
ability of particles with the wake-mediated interactions to 
form active units, the unusual melting and the unique dif¬ 
fusive behavior make such systems interesting for many 
fields of research. We encourage scientists in the field 
of colloidal dispersions or complex plasmas to realize the 
experiments where these theoretical predictions can be 
verified. 
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